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Abstract. Results are presented from general relativistic numerical computa- 
tions of primordial black-hole formation during the radiation-dominated era of the 
universe. Growing-mode perturbations are specified within the linear regime and 
their subsequent evolution is followed as they become nonlinear. We use a spher- 
ically symmetric Lagrangian code and study both super-critical perturbations, 
which go on to produce black holes, and sub-critical perturbations, for which the 
overdensity eventually disperses into the background medium. For super-critical 
perturbations, we confirm the results of previous work concerning scaling-laws but 
note that the threshold amplitude for a perturbation to lead to black-hole forma- 
tion is substantially reduced when the initial conditions are taken to represent 
purely growing modes. For sub-critical cases, where an initial collapse is followed 
by a subsequent re-expansion, strong compressions and rarefactions are seen for 
perturbation amplitudes near to the threshold. We have also investigated the 
effect of including a significant component of vacuum energy and have calculated 
the resulting changes in the threshold and in the slope of the scaling law. 
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1. Introduction 

Cosniological structure formation is tliouglit to liave resulted from tlie growtli and 
evolution of small perturbations initiated at the time of inflation (see 1 and references 
therein). Inflationary models give rise to a spectrum of fluctuations on scales larger 
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than the cosmological horizon which then start to re-enter the horizon in the radiation- 
dominated era. Primordial black holes (PBHs) could be formed at this stage in extreme 
cases where the fluctuation amplitude exceeds a critical threshold value |21 El ^ . The 
masses of these PBHs could, in principle, span many orders of magnitude, from the 
Planck mass up to the horizon mass at the time of equivalence between radiation and 
matter (~ 10^ years after the Big Bang). 

Sufficiently small PBHs would have emitted a significant amount of Hawking 
radiation and, according to the standard picture, would have evaporated completely 
away before now if their mass was less than ^ 10^^ g. It has been suggested that 
extremely small PBHs of ^ 10^ — 10^ g, formed in the very early Universe, might have 
been the cause of baryogenesis while PBHs of 10^^ g evaporating now might 
explain a class of observed very short-duration gamma ray bursts |S1 17| |Hj and some 
cosmic rays |51 ^]. It has also been suggested that the evaporation of small black 
holes might not continue all the way to zero mass but might stop in the region of 
the Planck mass, and that surviving PBH remnants might explain some or all of the 
current cold dark matter ^2 ^1 E| . 

Larger PBHs with masses well above ^ 10^^ g could be observable by means of 
gravitational effects such as their contribution to the cosmological density parameter 
| 14l [T5| , gravitational radiation emitted from coalescing binary PBH systems ^| or 
micro-lensing (PBHs were suggested as possible candidates for some of the objects 
observed by the MACHO Collaboration ^2] and subsequent surveys, although the 
idea that they might comprise a significant fraction of the halo dark matter is now 
no longer favoured |18[ 119 1). All of these aspects have been studied and constraints 
obtained both on the fraction of matter in the Universe that could now be in the form of 
PBHs and also on the spectral index of fluctuations on small scales |2()l 1211 1^ 1281 124j . 

In this article, we present results from general relativistic numerical simulations 
of PBH formation in the background of an expanding universe in the radiation- 
dominated era. There have been a number of previous studies of this type and our 
aim here is to re-visit the subject area, highlighting some features which we think 
are important and interesting. A particular aspect of our study is that we use initial 
perturbations representing growing-mode overdensities with length-scales larger than 
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the horizon and stiU within the hnear regime. Their evolution is then followed as they 
subsequently become nonlinear. A convenient parameter for measuring perturbation 
amplitudes is the fractional mass-excess within the overdense region (denoted by 6); 
if this is greater than a critical value Sc (super-critical perturbations), the following 
nonlinear evolution will result in formation of a black hole, whereas if it is smaller 
than 6c (sub-critical perturbations), the overdensity will eventually disperse back 
into the surrounding medium. Determining the value of Sc is very important for 
the cosmological considerations mentioned above. 

Following the earliest papers on this subject |21 El 0j Carr (1975) [Tl) 
carried out a quantitative study based on a simplified model consisting of an 
overdense collapsing region, described by a closed Friedmann-Robertson- Walker 
(FRW) spacetime, surrounded by a spatially-flat FRW expanding background. For 
a radiation-dominated universe, his calculation gave Sc = 1/3 and the black holes 
formed had masses Mbh of the same order as the horizon mass at the time when the 
fluctuation first entered the horizon. Nadezhin, Novikov & Polnarev (1978) |2S] carried 
out the first detailed numerical study of PBH formation using a hydrodynamical 
computer code similar to those of May & White (1966) and Podurets (1964) 
| 27| using a "cosmic-time" coordinate with a diagonal metric which reduces to a form 
similar to that of the FRW metric in the absence of perturbations. A well-known 
difficulty with this approach is that in a continuing collapse, singularities typically 
appear rather quickly and stop the computation before the black hole formation is 
complete. In pSj , this difficulty was overcome by using an early form of excision with 
the evolution being stopped in the region where the singularity would appear. The 
qualitative features of the earlier picture were basically confirmed but it was found 
that the PBH masses were always much smaller than the horizon mass. 

Shortly afterwards, Bicknell & Henriksen (1979) |2S| carried out related 
calculations using a method based on integration along hydrodynamical characteristics 
which avoids the problems associated with the appearance of singularities. They used 
rather different initial data from that of 25 and found formation of black holes with 
masses of the same order as the horizon mass (or greater in cases where the overdensity 
in the initial perturbation was not compensated by a surrounding under-dense region) . 
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They noted the appearance of both ingoing and outgoing compression waves during 
evolutions leading to black hole formation. 

More recently, Niemeyer & Jedamzik [2211201 made further numerical calculations, 
particularly focusing on the relevance of scaling-laws for PBH formation. They found 
that Mbh follows a power law in {5 — 5c) when the latter is sufficiently small, which 
is a similar behaviour to that seen in critical collapse by Choptuik (31| and many 
subsequent authors (see Evans & Coleman [31] and the review by Gundlach (32] )■ They 
started from initial perturbations specified at the moment when the overdensity enters 
the horizon and then computed the subsequent evolution. They used a null-slicing 
code, following the formulation by Hernandez & Misner (1966) 34 , and obtained 
6c — 0.7 for each of the three types of perturbation profile which they studied. 

While Niemeyer & Jedamzik [30] demonstrated the existence of a scaling law for 
PBH masses down to around one tenth of the horizon mass, they did not investigate 
smaller masses and so it was not possible to determine whether the scaling law was 
likely to continue down to vanishingiy small masses (when 5 — )■ (5c) as in type II critical 
collapse. In fact, the calculations become very challenging from a numerical point of 
view when 5 is very close to 5c because of the appearance of strong shocks and deep 
voids outside the region where the PBH is forming. Hawke & Stewart (2002) 
addressed this problem using a sophisticated purpose-built code which allowed them 
to make calculations for values of 5 closer to 5c and to handle the strong shocks which 
occur in these cases. They found that the scaling law does not continue down to very 
small values of {5 — 5c) but rather reaches a minimum value for Mbh around lO^** of 
the horizon mass, with the limit resulting from the behaviour of the shocks produced 
in nearly critical collapse. 

Shibata & Sasaki (1999) PZI presented an alternative formalism for studying 
PBH formation using constant mean curvature time slicing and focusing on metric 
perturbations rather than density perturbations. They emphasised the importance 
of using initial data which can be directly related to perturbations arising from 
inflation. Their formulation was not restricted to spherical symmetry (as had been 
the case for the previous authors mentioned here) but they presented results only 
from spherical calculations. They located the threshold perturbation amplitude for 
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PBH formation (in terms of the metric perturbation) and found that this varied 
considerably depending on the density of the medium surrounding the density peak. 
They concluded from this that it is probably important to take into account spatial 
correlations of density fluctuations when considering PBH formation. 

It is not straightforward to make the link between results from this type of 
calculation and ones from the more standard approach focused on density fluctuations. 
However, this was addressed in a recent paper by Green et al. (2004) They 
calculated the PBH abundance produced from two different fluctuation spectra, using 
peaks theory together with the threshold criterion of Shibata & Sasaki |37] . They then 
compared the results of this with ones obtained from a standard calculation based on 
a Press-Schechter-like approach and using density perturbations. They found that the 
Shibata & Sasaki results are consistent with ones using a density perturbation if Sc 
lies in the range 0.3 < Sc < 0.5 and they noted the discrepancy with the results of |30j . 

As mentioned previously, our aim in the work described in this paper was to 
re-visit some issues arising from these earlier calculations. In particular, we wanted 
to check on the scaling laws found in [301 and to investigate the effect on them both 
of starting the initial perturbations earlier, as pure growing modes within the linear 
regime, and also of including vacuum energy (as a cosmological constant term). For 
our calculations, we used a null-slicing code, modifled from that of Miller & Motta 
(1989) [321, which implements an essentially identical method to that used in [30] 
(but with one key technical difference which we discuss in Section 2). This is not as 
sophisticated in its treatment of the hydrodynamics as the Hawke & Stewart code 

but is satisfactory for our purposes and has some useful advantages particularly 
as regards the specification of initial conditions. We have examined both super- 
critical and sub-critical perturbations. The most notable results are: (i) we find good 
agreement with the results of |30| when we use the same initial data; (ii) we obtain 
Sc — 0.45 when we use initial data specified as purely growing-mode perturbations, 
consistent with the results of Green et al. (Hi) when a cosmological constant 

term A is included, Sc increases linearly with increasing A; (iv) subcritical collapse 
with S sufficiently close to Sc produces strong compressions and rarefactions before 
the overdensity subsides back into the surrounding medium. 
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The organisation of the paper is as follows. In Section 2, we review the 
mathematical formulation of the problem and the calculation method used; Section 
3 presents results from our calculations of PBH formation and sub-critical collapse; 
Section 4 contains discussion and conclusions. We use units for which c = G = 1 
except where otherwise stated. 

2. Mathematical formulation & calculation method 

For our calculations we have followed a procedure very similar to that of Niemeyer & 
Jedamzik PU]. We will not repeat their full discussion of the method but nevertheless 
it can be useful to make some further comments about it here. 

As with most of the other literature on this subject, we are restricting attention 
to spherical symmetry, which very greatly simplifies the calculations, and we have 
used the formulations of the relativistic hydrodynamical equations given by Misner 
& Sharp (1964) 40^ and Hernandez & Misner (1966) Both of these are 

Lagrangian formulations, the first using a diagonal metric (with the time referred to 
as "cosmic time" which reduces to the familiar FRW time coordinate in the absence of 
perturbations) , and the second using an outward null slicing where the time coordinate 
is an "observer time" (the clock time as measured by a distant fundamental observer). 
The cosmic-time formulation is particularly simple and has the advantage of using 
a slicing which many people find intuitive; this was the approach used by May & 
White (1966) | l26i in their classic paper studying spherically-symmetric gravitational 
collapse. However, as mentioned above, this approach has a well-known drawback for 
studying black hole formation in that singularities are typically formed rather early in 
calculations of continuing collapse and it is not then possible to follow the subsequent 
evolution. The outward null slicing approach is particularly convenient for calculations 
involving black hole formation in spherical symmetry: anything which could not be 
seen by a distant observer (e.g. singularity formation) does not occur within the 
coordinate timespan, while all observable behaviour can be calculated. This is, in some 
sense, the optimal approach for studying black hole formation in spherical symmetry as 
seen by an outside observer (being linked directly to potential observations) although 
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more sophisticated slicing conditions have advantages for calculations away from 
spherical symmetry. 

Following the introduction of the observer-time approach by Hernandez & Misner 
in 1966 . it was implemented soon afterwards in unpublished calculations. A brief 
presentation of some results was given by Miller & Sciama (1980) 41, and a full 
discussion of the technique used and of results obtained was given subsequently by 
Miller & Motta (1989) A problem with the use of this method concerns the 

satisfactory specification of initial conditions which is not natural to do on a null slice. 
Because of this, Miller & Motta made a preliminary calculation using the Misner 
& Sharp 40 formulation in order to construct data on an outgoing null slice from 
initial conditions specified on a space-like slice; the observer-time calculation then 
proceeded from the null-slice data constructed in this way. Subsequently, Baumgarte, 
Shapiro & Teukolsky (1995) 021 made calculations using a similar technique and it 
was these which were used as a reference point by Niemeyer & Jedamzik |30| . 

For the calculations of PBH formation reported in this paper, we have used a 
modification of the Miller & Motta code [3^1 (designed for studying the collapse of 
an isolated object surrounded by vacuum) together with elements from the codes by 
Miller & Pantano (1990) ^ and Miller & RezzoUa (1995) gH (designed for following 
phase-transition bubble growth within an expanding universe). 

In the remainder of this section, we give a brief overview of the equations used 
and of our numerical code. The reader is referred to the papers quoted above for 
further details. 

2.1. The Misner- Sharp equations 

For calculations in spherical symmetry, it is convenient to divide the collapsing matter 
into a system of concentric spherical shells and to label each shell with a Lagrangian 
co-moving radial coordinate which we will denote with r. The metric can then be 
written in the form 

ds^ = -a^ dt^ + dr^ + {d9^ + sin^ 0dip^) , (1) 
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where R (the Schwarzschild circumference coordinate) , a and b are functions of r and 
the time coordinate t. This was the form used Misner & Sharp (1966) |4U) . 

For a classical fluid, composed of particles with nonzero rest-mass, it is convenient 
to use the rest-mass fi contained interior to the surface of a shell (or, equivalently, the 
baryon number) as its co-moving coordinate r. For the case of a radiation fluid (as 
studied in this paper), rest-mass and baryon number are not available as conserved 
quantities to be used in this way but a similar procedure can still be followed by 
introducing the concept of a conserved number of unit co-moving fluid elements (Miller 
&: Pantano 1990) Denoting a "relative compression factor" for these fluid elements 
by p (equivalent to the rest-mass density in the standard treatment), one then has 

dfi^ AirpR^bdr, (2) 



and identifying fj, and r then gives 
1 



47ri?2p- 

Following the notation of "lUj, we write the equations in terms of the operators 



(3) 



D.^Ui), (4) 



a 



dt 



and applying these to R gives 

DtR = U, (6) 

DrR = r, (7) 

where U is the radial component of the four-velocity in the associated Eulerian frame, 
using R as the radial coordinate, and F is a generalisation of the Lorentz factor. 

We are mostly dealing with processes occurring in the radiation dominated era 
of the Universe for which the equation of state of the matter can be written as 

P=le, (8) 

where p is the pressure and e is the energy density. For one-parameter equations of 
state of the form p ~ p{e), the system of Einstein and hydrodynamic equations can 
then be written as: 



DtU 



F M 



(9) 



Computations of primordial black hole formation 9 
Dtp^-^Dr{R^U), (10) 

Ae - ^Ap, (11) 
P 

Dra = —DrP, (12) 

e + p 

DrM = inVeR^, (13) 

where M is a measure of the mass-energy contained inside radius ^ and F can be 
calculated either from Q or from the constraint equation 

r2 = i + t/2_^. (14) 

R 

2.2. The Hernandez- Misner Equations 

Because of the problems mentioned earlier, Hernandez & Misner IM' introduced the 
concept of "observer time" , using as the time coordinate the time at which an outgoing 
radial light ray emanating from an event reaches a distant observer|. In the original 
formulation, this observer was placed at future null infinity but for calculations in an 
expanding cosmological background we use an FRW fundamental observer sufficiently 
far from the perturbed region to be unaffected by the perturbation. Along an outgoing 
radial null ray wc have 

adt = bdr, (15) 

and we define the observer time u by 

f du = adt — bdr, (16) 

with / being an integrating factor which needs to be determined. In terms of this, the 
metric becomes 

ds^ = -/2 du"^ - 2/6 dr du + R^ {d0^ + sin^ ddip"^) , (17) 
which is no longer diagonal. The operators equivalent to I0J and (jSJ are now 

I We note that a somewhat similar approach, but based on a double null foliation, has recently been 
used by Harada et al. BBi for studying some different aspects of PBIf formation. 
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where D^. is the radial derivative across the null slice and the corresponding derivative 
across the Misner-Sharp space-like slice is given by 

Dr = Dk - Dt. (20) 

The observer-time equations replacing the cosmic-time ones - (|13f) are then: 

1 



DtU = - 



1 



r M 

J^)D,p+-+inRp 



DM 



2UT\ 



DtP^ 
Dte^ 
Dkf^ 



DtU^DkU 

e + p^ 



2ur 



R 



DkU + 



Dtp, 
M 



i?2 



47ri?p , 



DkM ^ A-KR^[eT ~pU], 



(21) 
(22) 

(23) 

(24) 
(25) 



where Cg — \J (dpjde) is the sound speed, which is equal to l/\/3 in the present case. 
The quantity F is given by equation (|14|l . as before, and also by 



r = DkR-U. 



(26) 



Using equations (24), (25) and (14), it is possible to derive the following alternative 
equation for /: 

-(T + uy 



D, 



f 



= -A-KR{e+p)f . (27) 

In calculations concerning collapse of an isolated object surrounded by vacuum in an 
asymptotically-flat spacetime, the observer time is taken to be the clock time of a 
static observer at future null inflnity and so (F + U)/ f = 1 at the location of that 
observer (since T = \, U — Q and / = 1 there). It then follows from equation (|27|l 
that (F + U)/ f = 1 also at the surface of the collapsing object, since the right hand 
side of (|27|) is zero in vacuum. The condition 



/ - r + {/, 



(28) 



at the surface is used as a boundary condition for / and the values of / internal to 
that are then calculated from equation (|24|l . 
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The general condition for a trapped surface is D^R < 0. With outgoing null 
slicing, ZJfci? should be reached only asymptotically in the future, accompanied 
by the lapse / going to zero, and DkR should never become negative. In practice, care 
is required in order to achieve the exact synchronization of DkR — > and / — > in 
a numerical solution where the equations are discretized; if the synchronization is not 
achieved, negative values of DkR do appear and the evolution becomes unphysical. 

In the case of an isolated collapsing object surrounded by vacuum, using boundary 
condition (|28|l together with equation H24(l ensures the correct behaviour. However, 
for the present situation, where the surroundings are not vacuum and the spacetime 
is not asymptotically flat, it is necessary to proceed in a different way. We are 
wanting to synchronise the "observer time" with the clock time of a co-moving FRW 
fundamental observer at the outer edge of the grid (setting / = 1 there) and then to 
calculate / elsewhere using this as a boundary condition. For doing this, we found 
it essential to use equation (|27|1 . which guarantees synchronisation of DkR and 
/ — > 0, rather than equation (|24|l . which always eventually gave rise to unphysical 
behaviour with DkR < 0. We think that this is a crucial point in using an outward 
null slicing technique for any situation regarding black hole formation within non- 
vacuum surroundings. It will apply equally to calculations of black hole formation 
from core-collapse of high-mass stars. 

2.3. The calculation method 

As mentioned above, our calculations of PBH formation have been made using an 
explicit Lagrangian hydrodynamics code based on that of Miller & Motta (1989) 
but with the grid organised in a way similar to that in the code of Miller & 
Rezzolla (1995) 01], which was designed for calculations in an expanding cosmological 
background. The reader is referred to those papers for full details of the methods 
used and we repeat just some main points here. The method described by Niemeyer 
& Jedamzik 30 , following from the Baumgarte et al. code is similar in most 
respects. In order to achieve good grid coverage, we have used a composite prescription 
for the grid spacing, with A/i increasing exponentially going outwards through the 
inner and outer parts of the grid but remaining constant in the intermediate region. 
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We start our perturbations in the linear regime with length-scales larger than the 
cosmological horizon radius Rh, and have the grid reaching out to 10 Rh with around 
2000 grid points. 

For collapses leading to black hole formation, our calculations proceed in two 
stages: first, initial data is specified on a space-like slice at constant cosmic time, 
specifying the energy density e and the four-velocity component U as functions of R 
at an initial time ti. This data is then evolved using the Misner-Sharp equations of 
Section 2.1 so as to generate a second set of initial data on a null slice (at constant 
observer time). To do this, an outgoing radial light ray is traced out from the centre 
and parameter values are noted as it passes the boundary of each grid zone. The second 
set of initial data, constructed in this way, is then evolved using the Hernandez- Misner 
equations of Section 2.2. For sub-critical cases, which end eventually with dispersal of 
the perturbation into the surrounding uniform medium, we have continued with the 
Misner-Sharp cosmic-time approach throughout. 

The unperturbed background model is taken as a spatially-flat FRW model, for 
which r = 1 (giving U = ■\/2M/R) with e(r) = constant at any particular cosmic 
time t. The perturbations of e and U are then superimposed on this background in 
the way described in the next section. During the following evolutions, the metric 
functions a (for cosmic time) and / (for observer time) are set equal to unity at the 
outer boimdary of the grid, thus synchronising the cosmic and observer times with 
local clock time there as measured by the local co-moving FRW fundamental observer. 

3. Description of the calculations 

In this section, we describe our calculations for both super-critical and sub-critical 
perturbations. 

3.1. Initial Conditions 

As initial conditions for our calculations, we use spherically symmetric density 
perturbations superimposed on a uniform FRW background. The perturbations are 
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specified in terms of the dimensionless quantity 

6e{R) = (29) 

where e& is the uniform background density given by eb = 3/327rt at any particular 
time t. Niemeyer & Jedamzik |3U) used three different types of expression for 5,, (shown 
in Figure 1 of their paper): 



Gaussian: 



Mexican hat: 



Polynomial: 



64R) = Acxp[-j^^], (30) 



^e(i^)^A(l-^)exp(-^), (31) 



6e{R) = { (32) 

if i? > V^ILRh) 

where the parameters A and L set the amplitude and length-scale of the perturbation. 
For the Mexican- hat and polynomial perturbations, the excess energy in the overdense 
region is exactly balanced by the deficit in the outer underdense region (i.e. 

AnSeR^dR — 0), whereas the Gaussian ones have only an excess, decreasing 
asymptotically to the background value. The latter is not very satisfactory for 
cosmological perturbations and so we concentrate here on the first two types. 

We follow the previous literature in also using an "integrated" perturbation 5 
which represents the mass-energy excess in the overdense region with respect to that 
in a corresponding uniform solution: 

n-Ro 

AirSeR'^dR 



S^'-^ , (33) 

where Ro is the radius of the overdensity. (In the case of the Gaussian, Ro is defined 
as the radius at which has fallen to 1/e^ of its value at i? = 0.) 

We are mostly starting with perturbations which are still on length-scales much 
larger than the horizon scale [Ro/ Rh > 5] and are well within the linear regime 
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[typically 6 ^ 10^^]. Setting initial conditions within the linear regime makes it easy 
to specify consistent density and velocity perturbations representing a purely growing 
mode: Se grows linearly with time and the associated velocity perturbation is given 
roughly by 

Su{R)=-j, (34) 

(coming from formulae in ,45,), where 6u is defined in the same way as for 6e in 
equation (f^ : 

S^iR) ^ (35) 

with Ub{R) = HR being the velocity field of the background Hubble flow. (We 
are here considering adiabatic perturbations with a Newtonian gauge, in which 
the perturbations do not generate off-diagonal components of the metric.) The 
corresponding decaying modes have decreasing amplitudes and, if excited by some 
process, would disappear again rather rapidly. Of the perturbations produced from 
inflation, only growing modes are of interest at the times which we are considering 
here. 

3. 2. Evolution of supercritical perturbations 

In this sub-section, we describe results from representative evolutions leading to black 
hole formation, carried out using the null-slicing formulation. 

Before moving on to new calculations, we first needed to check on whether our 
code did reproduce the results of Niemeyer & Jedamzik |30| when we used their 
choice of initial conditions and a simple exponential grid similar to theirs. We found 
extremely close agreement. For perturbations with the initial S only slightly larger 
than the critical value Sc, the masses of the black holes produced, Mbh, follow a 
scaling law: 

Mbh = K{6 - (t^ ), (36) 

where K and 7 are constants and MuitH) is the cosmological horizon mass at the 
horizon-crossing time Ih (i.e. when Roitn) = RnitH))- The type of behaviour given 
by H3t)|l is familiar from the literature on critical collapse (which, however, is generally 
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considering collapse under simpler circumstances and not within the context of an 
expanding universe). Our results match very closely with those of fSTO j , with ~ 0.67 
for the Mexican-hat profile and ~ 0.71 for the polynomial and Gaussian profiles and 
with 7 between 0.36 — 0.37 in each case. These values for 7 are very close to the value 
0.356... calculated semi-analytically for the same equation of state within the standard 
critical collapse scenario |46) . 
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Figure 1. Scaling behaviour for Mbh as a function of (5 — Sc) calculated 
for growing-mode Mexican-hat perturbations specified within the linear 
regime. The filled circles refer to the standard calculation discussed in 
section 3.2, while the open circles are for a calculation including a non-zero 
cosmological constant A, as discussed in section 3.3, giving y = 3.0 x 10~^. 

Following these initial calculations to reproduce previous results, we then carried 
out furtlier ones starting from growing-mode perturbations specified within the 
linear regime and with Icngtli-scales larger than the cosmological horizon. These 
perturbations were evolved with the Misner-Sharp code until the moment when they 
entered the horizon and the current value of 5 was then calculated for use as our 
measure of perturbation amplitude in the discussion of scaling laws. Having done 
this, we switched to the Hernandez-Misner code for completing the calculation. 
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Plotting the eventual black hole mass against {S — Sc), we obtained scaling 
curves very similar to those of 30 (and with almost identical values of 7) but with 
substantially different values for 6c'. for Mexican-hat perturbations, we found Sc — 0.43 
and for polynomial perturbations, 6c — 0.47. Our scaling-law results for Mexican-hat 
perturbations are shown as the A = curve in Figure ^ The reason for the changed 
values of 6c is clear: in |30) the initial perturbations, specified at the horizon-crossing 
time, had part of their amplitude contributed by a decaying-mode component which 
then rapidly decreased leaving only the growing-mode part visible. Only the part 
of the perturbation amplitude corresponding to the growing mode is relevant for the 
black hole formation and so the effective 6 is smaller than that calculated in fSU) . 

Noting the work of Hawke & Stewart (2002) and the fact that the horizon- 
scale is an important length-scale in the problem, we do not expect that the linear 
scaling law will continue to indefinitely small values of Mbh and {S — Sc) but, instead, 
that it will level off at some minimum value of Mbh- Confirming this behaviour in the 
case of growing-mode perturbations starting in the linear regime is of great interest 
but would require refinement of our present code to increase resolution and improve 
the capacity for handling strong shocks. 

Figure |21 shows some more detailed results for a particular representative case 
within the linear scaling regime. This run starts from a growing-mode Mexican-hat 
perturbation with Ro/Rh = 5 giving (6 - 6c) = 2.37 X 10"^ at horizon crossmg, 
which leads to formation of a black hole with Mbh ~ OAAld Mnitn)- The top 
two panels show the evolution of the lapse / and the corresponding behaviour of 
the fluid worldlines. The interpretation of the collapse of the lapse is particularly 
clear when an observer-time formulation is used: as / 0, the redshift of outgoing 
signals increases and the evolution as seen by a distant observer becomes frozen, 
corresponding to black hole formation (see also the small inset where log / is plotted) . 
Note, however, that strictly "black hole formation" occurs only asymptotically in the 
future according to this formulation. In the plot for the worldlines, one can see the 
separation between the matter which goes to form the black hole and the matter which 
continues to expand, with a semi-evacuated region being formed between them. Note 
that some of the outer material first decelerates but then accelerates again before 
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Figure 2. A typical evolution leading to black hole formation: the initial 

perturbation had a Mexican-hat profile and gave {S - (5c) = 2.37 X 10"^ at 
the horizon crossing time. The top left-hand panel shows the behaviour of 
the lapse function (the time sequence of the curves goes from bottom to top 
on the right hand side); the top right-hand panel shows the fluid-clomont 
worldlines (the time is measured in units of the horizon-crossing time tn)- 
The bottom left-hand panel shows the profile of 2M/R at different times, 
with the inset showing the approach of the maxinmin value of 2M / R 1; 
the bottom right-hand panel shows the corresponding evolution of the mass- 
energy (in both of these panels, the time sequence of the curves goes from 
top to bottom on the right hand side). 
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crossing this semi-evacuated region to fall onto the black hole. The bottom left-hand 
panel of Fig. |21 shows the behaviour of the ratio 2M/R, plotted against R at successive 
times. The event horizon corresponds to the asymptotic location of the outermost 
trapped surface, where DkR — and R — 2M. For the present purposes, we need 
an operational definition for calculating Mbh, bearing in mind that the black hole is 
only formed asymptotically and that further material may continue to accrete. The 
inset in the bottom left-hand panel of Fig. |2] shows the approach of the maximum 
value of 2M/R 1 : our operational definition for Mbh is to set it equal to the value 
of M at the maximum of 2M / R when (1 — 2M/R) there first becomes smaller than 
10^*. The bottom right-hand panel shows a corresponding plot for M against R. It 
can be seen that the profiles for M become very flat just outside the black hole region 
at late times, a consequence of the very low densities being reached there (less than 
10^* of the background density at the horizon-crossing time). The small inset shows 
the continuation of this figure up to larger radial scales. For calculations with S closer 
to 6c, the rarefactions formed become increasingly deep, and one sees strong shock 
waves appearing at the outer edge of the under-dense region. 

3.3. Evolution of super- critical perturbations when A > 

We were also interested to investigate the effect for PBH formation of including a 
cosmological constant large enough to affect the dynamics. We recognise that this 
is a highly idealised scenario since the present-day cosmological constant would have 
had negligible effect in the early universe and other vacuum energies present after 
inflation are unlikely to have been constant in time (e.g. quintessence). However, it 
is of conceptual interest to find out what the behaviour would be in this hypothetical 
case. A cosmological constant A is equivalent to a false vacuum with energy density 
Cv — A/Stt and pressure py — —A/Stt. Its effect can be included by adding these terms 
onto the standard energy density e and pressure p wherever those appear. 

A positive A eventually causes the expansion of the universe as a whole to start 
accelerating and acts against the growth of overdensities, while a negative A would aid 
general collapse. (The equations governing the background expansion when A 7^ are 
summarised in the Appendix.) As the background density decreases with time, the 
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cosmological constant becomes progressively more important until, when becomes 
greater than e, the deceleration is reversed and becomes an acceleration [see equation 
HA.2|) ]. It is convenient to introduce a quantity y, the ratio between the vacuum energy 
and the total energy in the uniform background 

y^^, (37) 

e + Cy 

which can then be written as 

y^^AMfj, (38) 

since Mh = ^T:RH^{e + e„) with Rh — 2Mh (Note that this type of relation holds 
for the cosmological horizon in the same way as for a black hole event horizon) §. In 
the following, we will use y as a general measure of the importance of the A term, 
with Mh being measured at the horizon-crossing time for a perturbation with S — dc- 
The influence of A during formation of a black hole of mass Mbh can similarly be 
characterised by the quantity AM^jj and so, for a given A, is greatest for large black 
holes. 

In making computations with A > 0, it is particularly important to start the 
calculation at an early time when the perturbation has a length-scale larger than the 
horizon. For appreciating its effects, one wants A to be sufficiently large so as to make 
a significant difference for the collapse, but not so large that it creates a problem for 
constructing the null-slice initial data. By starting the calculation sufficiently early, 
this can be achieved although the values of y which we will be considering are all very 
small (in the range 10^^ — 10^^). 

The qualitative picture of collapses leading to black hole formation is not changed 
very greatly by the presence of a A term but there are significant differences in the 
parameters of the scaling law. We started all of these calculations with perturbations 
at five times the horizon scale (i.e. Rq/Rh — 5). 

The impact of A on the scaling law can be seen in Figure^] 7 decreases for A > 
and for sufficiently small A is found to follow a linear relationship 

7(A) -7(0)- 8.3 y. (39) 

§ Both are trapped surfaces. A black hole event horizon is the asymptotic location of the outermost 
trapped surface for outgoing light-rays whereas the cosmological horizon is the innermost trapped 
surface for incoming light rays. 
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The critical amplitude 6c increases with increasing A and also follows a linear 
relationship: 

6ciA)~Sc{0) + 0.98y. (40) 

This behaviour can be interpreted as follows: a positive A acts against collapse, so 
that corresponding black hole masses will be lower and the threshold amplitude Sc 
will be raised. For a given A, its influence is greater for larger black-hole masses than 
for smaller ones (cx M^^) and this gives rise to the observed decrease in 7. 

3.4- Evolution of subcritical perturbations 

For subcritical perturbations with d considerably less than 6c, the perturbation initially 
grows but then subsides back into the surrounding medium in an uneventful way. 
However, for perturbations with 6 sufflciently close to Sc, some very interesting 
behaviour is seen and we present results from a representative case of this in the present 
subsection. Our calculations for subcritical perturbations use only the Misner-Sharp 
code. 

The run presented starts with a Mexican-hat perturbation specified in the linear 
regime, with {5 — Sc = —3 x lO^'^). In Figure |2| the fluid worldlines are plotted and 
the main features of interest can already be seen from this. Figures 01 and [S] then show 
details of the evolution of the energy density e and the radial velocity U. Figure 01 
shows a sequence of snap-shots of these quantities, plotted as a functions of R, at key 
moments during the evolution. Figure El shows the time evolution of these quantities 
at three (comoving) locations: near the centre of the perturbation, at an intermediate 
region (mid-way through the collapsing matter) and at the edge of the grid where the 
fluid is unperturbed. 

The evolution can be summarized in terms of the following steps (see particularly 
Figures 13 and ^: 

(i) Initially, the perturbation has very small amplitude {Sc ^ 10~^) and its length- 
scale is five times the horizon scale; the perturbation amplitude then grows within 
the expanding fluid. The deceleration in the perturbed region is larger than that 
in the unperturbed region and its expansion lags progressively behind that of 
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Figure 3. Worldlines for a Mexican-hat perturbation with {S — 5c) = 
"3.0 X 10""^). This plot shows alternating collapse and expansion of the 
perturbed region while the outer material continues to expand uniformly. 
The "cosmic" time is measured in units of the time at horizon crossing. 

the outer matter until eventually it starts to re-contract shortly after horizon 
crossing. The maximum infall velocity reached is, however, rather small [see row 
(a) of Figure^ which is for a time considerably after horizon crossing, when the 
perturbation has become very nonlinear] . The infall can be clearly seen in Figures 
01 and 01 but is only just visible in the second frame of Figure El 

(ii) The contraction is not strong enough to produce a black hole and the fluid bounces 
out again [row (b)], expanding until it encounters the surrounding matter which 
did not participate in the contraction. A compression wave forms where the two 
regions of fluid meet, while the density becomes very low at the centre of the 
perturbation [row (c)]. 

(iii) The compression wave proceeds out into the surrounding material [row (d)] 
but also some matter is sent back into the middle of the rarefaction where it 




Computations of primordial black hole formation 22 




J I I I I I I I I I I I I I I I i_d I I I I I I I I I I I I I I I I I L 

5 10 15 5 10 15 



Figure 4. Plots of local quantities as functions of R/RnitH)'- the velocity 
U/c is shown in the left-hand column and the energy density e/eb in 
the right-hand column. The frames correspond to the following values 
of {t~to)/tH: (a) 7.02; (b) 25.92, (c) 31.67; (d) 33.64; (e) 40.11. Note 
that Rh is increasing with time and so points with R/RuitH) > 1 can be 
within the current horizon scale at times after horizon-crossing. 
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undergoes a second bounce which is much more extreme than the first with a 
very abrupt change of velocity in the central regions (as can be seen from Figure 
[SJ. Whereas the outward moving compression is damped geometrically as it 
proceeds to spherical surfaces with progressively larger areas, the inward-moving 
wave of material is geometrically amplified by the inverse process. The reason 
for the second collapse and bounce being more violent than the first is that while 
the first is a collapse of an ovcrdcnsity which is resisted throughout by internal 
pressure, the second is essentially the collapse of a "shell" with near vacuum inside 
it and is hence close to free-fall until just before the bounce. 

(iv) The compression wave formed by the second bounce propagates out into the 
surrounding medium following the first one [row (e)]. Both proceed to damp 
geometrically and eventually the medium returns to a uniform state. Note that 
the second compression wave is very steep fronted (see Figure 0J but is not quite 
a shock. It is likely that that genuine shocks would be seen for perturbations with 
S closer to Sc- 
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Figure 5. Evolution of the energy density e and radial velocity U at 
three (comoving) locations: near to the centre of the perturbation, at 
an intermediate region and at the edge of the grid where the fluid is 
unperturbed. Each quantity is measured in units of its initial value at 
the same comoving location. 
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It is useful to make some further comments about the plots in Figure [3 The 
left-hand plot shows the density normalized to its initial value at the same comoving 
location: this shows the variation of the local density against the background of the 
general decrease in density as the universe expands and allows one to see clearly the 
local contraction and expansion of the fluid. Initially, the perturbed region is also still 
expanding but it then starts to contract when the perturbation amplitude has reached 
a high enough value. The right-hand plot shows the local value of the radial velocity 
also normalized to its initial value at the same co-moving location. This has a very 
different behaviour from that of the energy density, with the inward velocity during 
the second collapse being very much greater than that during the first and the second 
bounce being far more abrupt. We see just two bounces in our calculation: after this 
the expansion of the universe prevents further ones, although we expect that more 
bounces would be seen for an initial perturbation with 6 closer to Sc- 

4. Conclusions 

We have presented results from calculations of PBH formation in an expanding 
universe during the radiation dominated era, following a broadly similar approach 
to that used by Niemeyer & Jedamzik 30 in previous investigations of this subject 
but with a number of key differences. Growing- mode perturbations were introduced 
into an otherwise uniform expanding medium and their subsequent behaviour was then 
followed. Those overdensities with amplitudes S greater than the critical value Sc give 
rise to black hole formation whereas smaller ones contract to a maximum compression 
and then bounce, eventually dispersing into the surrounding medium. Defining S to 
be the fractional mass excess within the overdensity at the time of horizon crossing, 
we find that 5c ranges between 0.43 and 0.47, depending on the perturbation shape. 
These values are consistent with the results of Shibata & Sasaki 37 , obtained with 
a different method, resolving a previously noted discrepancy .'iS . We have verified 
the results of (SO] concerning the scaling law relation between the black-hole mass 
Mbh and {S — Sc) when the latter is sufficiently small. The power of the scaling 
law 7 is essentially unchanged for our modified initial conditions, in contrast with 
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the large changes seen for 5c- We note, however, the important resuhs of Hawke & 
Stewart (2002) ^] who found that, for a class of perturbations specified in a non- 
linear regime at sub-horizon scales, the scaling law does not continue to arbitrarily 
small Mbh but levels off at around 10^'* of the horizon mass due to interaction with 
the surrounding medium. Our code does not yet have the capability of approaching 
sufficiently close to Sc to reveal this effect; it would be important for future work to 
establish to what extent the Hawke & Stewart result is dependent on the form of 
the initial conditions. 

We have investigated the effect on the scaling law of introducing a cosmological 
constant term A and found that positive values of A give rise to lower values of 7 and 
higher values of (5c, each varying linearly with A when the latter is sufficiently small. 
This can be understood in terms of the effect of a cosmological constant in opposing 
collapse. 

We also studied sub-critical collapses (with 6 < Sc). When 6 is close to Sc, we 
find that this can be a surprisingly violent process. The initial collapse and bounce 
are rather mild, being moderated by the internal pressure, but after the subsequent 
reflection from the surrounding medium and re-collapse (with matter falling into near- 
vacuum) , the second bounce can be very violent with a rapid velocity change and the 
sending out of a second compression wave into the surrounding medium following 
the one produced by the first bounce. After this, the possibility of further bounces is 
stopped by the overall background expansion for the cases studied (although additional 
bounces are expected for 6 closer to 6c)- 

There are a number of topics concerning PBH formation which may be of 
great interest for cosmology and which we are intending to study further, including: 
determination of the proportion of the matter in the universe going into black holes 
according to various scenarios for the perturbations coming from inflation (with the 
possibility of ruling out some scenarios); investigation of possible enhanced PBH 
formation at the time of phase transitions (see Jedamzik & Niemeyer 03); alternative 
possibilities for the formation of intermediate mass black holes or of the seeds for super- 
massive black holes at the centres of galaxies (as mentioned by Bicknell & Henriksen 
|2H1 hi 1979). Work is now in progress to investigate these. 
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Appendix A. Cosmological solution with non-zero A 

We present here some of the analytic expressions describing an expanding universe 
with A 7^ in the radiation-dominated era. These equations are certainly not being 
presented here for the first time but they are not easy to find in the usual cosmology 
textbooks and we think that it may be useful to present them together here. In this 
part only, we use physical units and do not set c = G = 1. 

First, we note the forms taken by the Friedman equation and the associated 

acceleration equation (we are taking the spatially flat case): 

a \ 8ttG Ac^ , . , 

, (A.l) 



a J 3c^ 3 
a\ SttG / Ac^ 



J 3c2 V SttG / ' 
where e is the energy density of radiation which scales as 



(A.2) 



e = e.(^) , (A.3) 

(here the subscript i refers to a fiducial initial time). Inserting this into ljA.l|l . we 
obtain the integral equation 

n(a,/a) —d{ — ] /-t 



SvrG Ac2 /a,\* -'o 



-ei 



3c2 ' 3 \a' 
If A > 0, the solution for the scale factor is 



ait) = ai ( 



SttGc, 



1/4 



smh I 2\l —ct 



Using l|A.5|) in IIA.l|) we get the Hubble parameter 



1/2 



(A.5) 



H{t)^^^ccoi\li^^j^cty (A.6) 

which, in the limit A ^ 0, reduces to the standard expression H{t) — l/2t. Inserting 
HA.6|I into IjA.ljl we get the expression for e{t) 



, , Ac^ 



sinh I 2\l —ct 



(A.7) 
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H 



hi 



(A. 



which we have used in the cosmic time code to calculate the initial time for the 
calculation. 

Finally, for completeness, we calculate the expression for the two cosmological 
horizons. From (|A.6|I we get a straightforward the expression for the Hubble horizon 
Rh = c/H, 



RH{t) 



tanh I 2\l —ct 



(A.9) 



For the particle horizon, the calculation is more complicated. From the definition 

cdt' 



Rh (t) = a (t) 



a{t') 



we get 



Rh (t) = 



sinh I 2\l —ct 



1/2 



du 



sinh ( 2\/ jcu 



1/2 



(A.IO) 



(A.ll) 



and, with the aid of integral tables, we get the final form of the solution 

-1 r / , — \ n 1/2 



Rh{t)= \ 2\j- 



sinh I 2\j —ct 



where F{(j)^ k) is an incomplete elliptic integral of the first type: 

d0 



Fi(b,k) 



with 



(l-Fsin^ 61)1/2' 



1 -sinh(2jfrf 



and 



arccos 



1 



1 + sinh 2 



fc = ^. 
V2 



(A.12) 



(A.13) 



(A.14) 



(A.15) 



When A < 0, relations (|A.5|I . HA.6|I and IjA.?!! are unchanged apart from replacing 
the hyperbolic functions by the corresponding trigonometric ones. This is consistent 
with the oscillating behaviour of a universe with A < 0, characterized by a sequence 
of expanding and contracting phases. 
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